Levy  searches  based  on  a  priori  information:  The  Biased  Levy  Walk. 
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Searching  for  objects  with  unknown  locations  based  on  random  walks  can  be  optimized  when 
the  walkers  obey  Levy  distributions  with  a  critical  exponent.  We  consider  the  problem  of  opti¬ 
mizing  statistical  searches  when  a  priori  information,  such  as  location  densities,  are  known.  We 
consider  both  spatially  dependent  exponents  and  biased  search  directions.  For  spatially  localized 
target  distributions  and  non-destructive  searches,  the  search  is  most  improved  by  biasing  the  search 
direction. 


The  best  statistical  strategy  for  efficient  searching  for 
randomly  located  objects  is  a  subject  of  ongoing  investi¬ 
gation  [1—4].  Levy  random  walks  are  known  to  outper¬ 
form  Brownian  (normal)  random  walks  [3]  when  the  pre¬ 
cise  location  of  the  targets  is  not  known  a  priori  but  their 
spatial  distribution  is  uniform.  Levy  flights  are  charac¬ 
terized  by  a  distribution  function 

n (I,)  ~  ija ,  (i) 

where  lj  is  the  flight  length  and  1  <  a  <  3  [2].  The 
lower  bound  of  1  is  required  in  order  to  have  a  normal¬ 
ized  probability  distribution  (1).  Such  statistical  search¬ 
ing  patterns  are  also  robust  as  they  result  in  space  fill¬ 
ing  paths.  For  exponent  values  exceeding  3,  the  central 
limit  theorem  implies  that  the  distributions  are  Gaus¬ 
sian.  The  optimally  efficient  searching  exponent  can  be 
computed  analytically  for  a  nondestructive  search  for  a 
spatially  uniform  target  distribution  [2].  In  the  limit  of 
large  mean  free  path  between  targets  one  obtains  a  =  2; 
for  destructive  searches  a  j  1  is  optimal  [3]. 

A  uniform  target  distribution  applies  to  problems  in 
biology  such  a  foraging  for  food  [5].  However,  for  hu¬ 
man  applications,  in  which  the  searchers  are  software, 
robots,  humans,  etc.,  there  is  often  presumptive  knowl¬ 
edge  and  there  are  many  approaches  to  best  use  this  in¬ 
formation.  For  example,  a  recent  study  of  U-boat  search 
and  destroy  missions  [6]  proposes  a  human-in-the-loop 
integrated  with  a  heuristic  optimization  method.  In  this 
paper,  we  propose  an  autonomous  search  method  based 
on  the  classical  Levy  method  but  applied  to  a  problem  in 
which  the  distribution  of  targets  is  spatially  non-uniform. 
We  consider  both  spatially  dependent  exponents  a  and 
biasing  the  choice  of  direction  for  each  successive  flight. 
For  simplicity  we  consider  a  Gaussian  target  distribution. 
Future  work  would  generalize  to  more  complex  target  dis¬ 
tributions. 

We  build  our  model  on  an  unbiased  Levy  search  with 
flight  lengths  satisfying  Eq.  (1),  in  which  the  searcher 
“walks”  the  path  as  opposed  to  making  an  instantaneous 
jump  [2].  In  this  model  the  searcher  behaves  as  follows: 
(a)  If  there  is  a  target  located  within  a  ‘direct  vision’ 
distance  rv,  then  the  searcher  detects  the  target,  (b) 
If  there  are  no  targets  within  a  distance  rv,  then  the 
searcher  chooses  a  direction  at  random  and  a  distance  lj 
from  a  Levy  distribution.  The  searcher  then  incremen¬ 


tally  moves  to  this  new  point,  constantly  searching  for  a 
target  site  within  a  distance  rv.  At  the  end  of  the  jump, 
the  searcher  picks  a  new  jump  distance  and  direction  and 
repeats  the  process.  Multiple  targets  could,  in  principle, 
be  detected  during  any  individual  jump  path. 

The  model  on  which  we  build  has  only  two  parameters, 
the  vision  distance  rv  and  the  exponent  a  of  the  Levy 
distribution.  The  vision  distance  is  a  fixed  parameter  for 
the  model  determined  by  the  physical  limitations  of  the 
searcher.  Although  we  develop  theory  and  algorithms  for 
an  infinite  domain,  for  practical  purposes  and  for  numer¬ 
ical  simulation,  we  restrict  to  a  finite  system  size.  In  this 
paper  we  contrast  two  methods  of  incorporating  a  pri¬ 
ori  information  into  the  standard  Levy  model:  the  first 
allows  the  distribution  exponent  a  to  depend  on  space. 
The  second  biases  the  direction  chosen  for  each  successive 
flight  path  (the  standard  model  has  equal  probability  to 
jump  in  any  direction). 

We  first  review  the  results  in  [2]  on  the  optimal  expo¬ 
nent  a  for  both  nondestructive  and  destructive  searches. 
Let  A  denote  the  mean  free  path  of  the  searcher  between 
successive  target.  The  mean  flight  distance  is 
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Here  we  assume  that  the  mean  distance  between  succes¬ 
sive  targets  is  a  fixed  A  which  thus  serves  as  a  maximum 
necessary  jump  length.  Likewise,  no  flights  are  less  than 
the  “vision  distance”  rv.  The  search  efficiency  function 
o(ct)  is  the  ratio  of  the  number  of  target  visited  to  the 
total  distance  traveled  by  the  searcher, 
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where  K  is  the  mean  number  of  flights  needed  to  travel 
between  two  successive  targets.  H  can  be  computed  as 

K  ~  (A/r„)Q-1for  destructive  searching,  (4) 

H  ~  (A/r„) for  nondestructive  searching,  (5) 


for  1  <  a  <  3.  In  the  case  where  targets  are  sparsely 
distributed,  A  rv,  then  substituting  Eq.  (4)  and  Eq. 
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(2)  into  (3),  one  finds  that  for  destructive  searching,  the 
mean  efficiency  has  no  maximum,  with  lower  values  of  a 
leading  to  more  efficient  searching.  Thus  for  destructive 
searching  the  optimal  choice  is  a  as  close  to  1  as  possible. 
For  nondestructive  searches,  replacing  Eq.  (4)  with  (5)  in 
Eq.  (3)  and  differentiating  with  respect  to  a ,  the  optimal 
efficiency  is  achieved  at 

a* (A,  rv)  =  argmaxa>17/(a,  A,r„)  =  2  -  <f(A,  r„),  (6) 

where  6  ~  l/[ln(A/r„)]2.  In  summary,  for  a  nondestruc¬ 
tive  search,  an  optimal  searching  exponent  a*  is  pre¬ 
dicted;  where  a*  is  a  function  of  the  mean  target  distance 
A  and  the  vision  radius  rv.  These  results  are  independent 
of  the  dimension  of  the  search  domain. 

Searching  in  one  dimension.  We  propose  two  meth¬ 
ods  for  constructing  a  spatially  dependent  exponent  a. 
The  first  method  directly  applies  the  above  arguments 
for  uniform  target  distribution  to  a  nonuniform  distribu¬ 
tion.  The  second  method  is  a  heuristic  model  motivated 
by  observational  data  from  simulation.  Both  models  re¬ 
quire  an  upper  bound  on  the  distance  between  targets  as 
a  function  of  position  in  the  spatial  domain.  Assume  a 
probability  distribution  of  target  locations  P(x)  ( P  >  0, 
f  P  =  1).  Then  the  expected  length  between  a  fixed 
point  x  in  the  search  domain  and  a  point  2  chosen  from 
the  target  distribution  is  given  by 

/OO 

\x  —  z\P(z)dz.  (7) 

-OO 

This  gives  an  effective  upper  bound  on  the  distance  to  the 
nearest  target.  One  strategy  for  improving  a  searching 
method  for  a  spatially  variable  target  distribution  is  to 
take  a  spatially  dependent  a  based  on  a  local  value  of  A. 
We  compare  two  contrasting  options  using  A  (a:)  above. 

Option  1A.  Use  a*  in  formula  (6),  plugging  in  for 
A  the  spatially  dependent  bound  computed  in  (7).  Note 
that  for  a  fixed  vision  radius  rv,  a*  increases  (to  the  value 
2)  as  A  increases.  A  larger  value  of  a  makes  the  search 
closer  to  Brownian  motion,  thus  keeping  the  searcher 
more  contained  in  a  region.  This  is  because  the  larger  the 
exponent  a,  the  smaller  the  average  jump  length,  with 
the  jumps  approximating  Brownian  motion  as  a  — >  3. 

Option  IB.  For  a  spatially  localized  probability  distri¬ 
bution,  we  might  heuristically  expect  that  a  better  search 
strategy  uses  a  smaller  a  where  targets  are  sparse  (larger 
A)  and  a  larger  a  where  targets  are  closer  together.  This 
way  a  searcher  would  take  smaller  hop  lengths  in  a  dense 
target  region  and  larger  hop  lengths  in  a  less  dense  re¬ 
gion. 

Using  the  mean  distance  to  a  target  computed  in  (7), 
we  propose  to  choose  a  according  to  a 2(2;)  = 

(l  +  =^kU)+3H(r„- A(*)).  (8) 

Here  H  denotes  the  Heaviside  function  which  is  one  for 
positive  arguments  and  otherwise  zero.  Note  that  Option 
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FIG.  1:  The  figure  contrasts  a*  (solid  line)  with  02  (dashed 
line)  for  a  Gaussian  distribution  P(x)  with  variance  1,  zero 
mean,  and  a  vision  radius  rv  =  0.5. 

IB  has  a  decreasing  as  A  increases.  Figure  1  contrasts 
a*  and  02  for  the  case  of  a  Gaussian  target  P  with  zero 
mean  and  unit  variance  and  a  vision  radius  rv  =  0.5.  For 
this  particular  choice  of  P,  A  can  be  computed  explicitly 

as  A(:r)  =  crerf  e-®2/2.  Note  that  a*  increases 

for  larger  x  while  a  2  decreases.  The  significance  of  this 
is  that,  for  the  Gaussian  distribution,  Option  1A  results 
in  shorter  jump  lengths  in  regions  of  lower  target  density 
while  Option  IB  results  in  longer  jump  lengths  in  these 
regions. 

Option  2.  We  propose  a  third  method  that  biases  the 
jump  direction  based  on  the  a  priori  information.  Af¬ 
ter  choosing  a  jump  length  from  the  Levy  distribution,  a 
choice  of  direction  is  made  using  information  about  the 
derivative  Px  of  the  a  priori  target  distribution.  Since 
we  wish  to  retain  the  space  filling  nature  of  the  search, 
we  choose  a  nonzero  probability  of  moving  in  either  di¬ 
rection.  We  propose  the  following  algorithm  in  which  (3 
denotes  a  number  between  0  and  ^  and  pr  and  pi  denote 
the  probabilities  of  jumping  to  the  right  and  to  the  left. 

•  Px  >  0,  (the  a  priori  distribution  is  increasing), 
then  choose  pr  =  1  —  /3,  pi  —  /3. 

•  Px  <  0,  (the  a  priori  distribution  is  decreasing), 
then  pr  =  /3  and  pi  =  1  —  j3. 

•  Px  =  0,  then  perform  unbiased  Levy  search. 

In  the  computations  below,  we  choose  (3  =  0.25.  In  one 
dimension,  this  method  is  very  similar  to  the  classical  bi¬ 
ased  random  walk  [7].  Here  we  use  local  spatial  informa¬ 
tion  to  determine  the  bias.  In  two  dimensions  (discussed 
below)  the  method  is  more  complicated. 

We  present  simulations  of  the  three  methods  above  for 
nondestructive  searching  on  a  fixed  box  of  length  200. 
The  vision  radius  rv  =  1.5  and  the  searcher  is  started 
randomly  in  the  domain.  Figure  2  shows  the  result  of 
each  method  as  a  function  of  N,  the  number  of  targets. 
Each  datapoint  is  an  average  over  1000  runs  of  100  jumps 
each.  Figure  2  shows  that  the  biased  Levy  search  out¬ 
performs  all  unbiased  searches.  The  variable  «2  (Option 
IB)  method  is  the  next  efficient  followed  by  the  standard 
unbiased  with  a  =  2.  Finally  the  simulations  that  used 
a  =  a*  (Option  1A)  are  least  efficient.  At  the  surface  this 
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FIG.  2:  One  dimensional  Levy  walks,  variable  a  compared 
with  biased  searching  for  a  Gaussian  target  distribution  with 
variance  1.  Efficiencies  as  a  function  of  the  number  of  targets 
N  averaged  over  1000  runs  of  100  jumps  each.  The  vision 
radius  rv  =  1.5.  The  symbols  denote  ■  =  a*,  A  =  a*  with 
bias,  ♦  =  unbiased  with  a  =  2,  •  =  02,  solid  line  =  biased 
with  a  =  2. 


is  surprising  as  a*  is  originally  derived  to  maximize  effi¬ 
ciency.  However  with  a  concentrated  target  distribution, 
this  choice  of  a  tends  to  have  a  larger  jump  length  which 
gives  a  higher  probability  of  “jumping  over”  the  concen¬ 
tration  of  targets  at  the  origin.  Figure  2  also  shows  that 
bias  search  outperforms  the  variable  a  approach.  We  ex¬ 
pect  that  higher  dimensionality  will  enhance  this  effect. 

Two  dimensions.  The  proposed  ID  search  methods 
can  be  extended  to  higher  dimensions.  Given  a  2D  prob¬ 
ability  distribution  P  of  target  sites,  the  expected  length 
to  a  target  can  be  computed  as 

/oo  /»oo 

/  \/{x  —  u )2  +  (y  —  v)2P(u,  v)dudv. 

-oo  J  —  oo 

(9) 

The  proposed  methods  in  Options  1AB  above  can  be 
directly  extended  two  dimensions  with  this  new  choice  of 
A. 

Extending  the  biased  Levy  search  requires  more  dis¬ 
cussion.  Given  an  a  priori  two  dimensional  target  dis¬ 
tribution  P(x,y),  we  choose  the  jump  direction  from  a 
continuum  as  opposed  to  left  or  right  as  in  the  one  di¬ 
mensional  case.  We  choose  an  angular  direction  from  a 
distribution  with  mean 


p  =  tan  1 


(  dyP{x,y)\ 
\dxP(x,  y) ) 


(10) 


In  the  case  that  VP(x,  y)  is  undefined  or  zero,  the  angle 
is  chosen  uniformly  in  [0,  2? r)  (i.e.  it  reduces  to  an  un¬ 
biased  Levy  step),  otherwise  we  choose  the  angle  from  a 
von  Mises  distribution  [8]  with  mean  p.  The  von  Mises 
distribution  for  points  on  a  circle  is  analogous  to  the  nor¬ 
mal  distribution  of  points  on  a  line.  Its  properties  are 


well  documented  in  [9].  An  angular  random  variable  9 
has  a  von  Mises  distribution  V M  (p,  n)  if  its  probability 
density  function  has  the  form 

$(fl)  =  -  }  ,  eKC°s(e-^,  -7 T<9<TT,  - TT<p<TT , 

2tt/0(k) 

where  Jo  (ft)  is  the  zeroth  order  modified  Bessel  function 
of  the  first  kind.  We  use  the  wrapped  Cauchy  method 
[10]  to  generate  variates  from  this  type  of  density.  The 
variables  p  and  k  are  analogous  to  the  mean  and  vari¬ 
ance,  respectively,  of  a  normal  distribution  on  the  line. 
We  choose  the  angle  of  direction  from  a  von  Mises  distri¬ 
bution  with  parameters  p  given  by  (10)  and  k  constant. 
If  the  gradient  of  P  is  zero  or  undefined,  then  9  is  chosen 
uniformly  in  [0,2n). 

Again,  we  simulate  and  compare  the  different  search 
methods.  As  before  we  choose  P  to  be  Gaussian  with 
zero  mean  and  unit  variance.  The  vision  radius,  rv  =  1.5, 
and  we  average  over  1000  runs,  each  truncated  after  100 
jumps.  A  smaller  domain  [—10,10]  x  [—10,10]  is  used; 
on  the  domain  [—100, 100]  x  [—100, 100],  the  majority  of 
searching  is  done  in  sparse  regions  preventing  meaning¬ 
ful  statistics  with  the  100  jump  cutoff.  On  the  smaller 
domain  [—10, 10]  x  [—10, 10]  we  can  achieve  a  similar  hit 
rate  for  targets  in  the  Gaussian  distribution  as  we  did  on 
the  [—100,100]  domain  in  one  dimension.  We  compute 
A (x,y)  numerically  on  the  lattice  with  spacing  A.,  and 
interpolate  as  needed  at  each  jump  point  in  the  simula¬ 
tion.  Figure  3  shows  the  efficiencies  as  a  function  of  tar¬ 
get  number  for  the  different  methods.  In  two  dimensions 
the  increase  in  efficiency  by  using  a  biased  search  (Option 
2)  is  more  pronounced  than  in  one  dimension.  All  three 
computations  without  bias  (fixed  a  =  2,  a  =  a.2  and 
a  =  <22)  performed  similarly.  The  biased  simulation  with 
a  =  2  outperforms  all  the  others,  including  the  biased 
simulation  with  a  =  a*. 


FIG.  3:  Two  dimensional  Levy  walks,  efficiencies  vs.  target 
number.  The  symbols  are  as  in  Figure  2,  ■  =  a*,  A  =  a* 
with  bias,  ♦  =  unbiased  with  a  =  2,  •  =  02,  solid  line  = 
biased  with  a  =  2. 
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In  two  dimensions,  Option  IB,  defined  by  Eq.  (8)  per¬ 
forms  less  well  than  in  one  dimension.  This  is  related  to 
the  fact  that  direction  is  more  important  than  in  the  one 
dimension.  This  effect  is  illustrated  in  the  trajectories 
seen  in  Fig.  4.  The  left  panel  denotes  the  trajectory  from 


FIG.  4:  The  left  panel  shows  a  biased  ( k  =  1)  two  dimensional 
trajectory.  The  right  panel  shows  a  dynamic  a  trajectory 
(using  a*  (Option  1  A) ) .  rv  =  0.5,  and  100  jumps  are  allowed 
before  cutoff.  The  step  size  is  1  and  every  5  steps  are  plotted. 

a  biased  (k  =  1  in  the  von  Mises  distribution)  Levy  flight, 
and  the  right  panel  the  trajectory  from  a  Levy  flight  with 
variable  a.  The  circles  are  the  level  sets  of  the  a  priori 
target  distribution.  When  the  trajectory  intersects  with 
the  boundary  of  the  search  domain  [—10, 10]  x  [—10, 10], 
the  jump  is  truncated  and  the  searcher  jumps  again.  The 
maximal  jump  length  is  restricted  to  be  the  minimum 
of  A(a;)  and  20.  With  these  constraints,  we  see  that  the 
majority  of  the  variable  a  search  is  spatially  near  the 
boundary  of  the  search  domain,  while  the  biased  search 
repeatedly  passes  through  the  peak  of  the  a  priori  target 
distribution.  The  variable  a  search  has  the  same  prob¬ 
ability  of  jumping  in  any  direction,  therefore,  when  the 
searcher  is  near  the  boundary,  there  is  a  50%  chance  that 
the  searcher  will  jump  towards  the  edge  of  the  search  do¬ 
main.  In  the  biased  case,  when  the  searcher  is  near  the 
boundary,  there  is  a  large  probability  that  the  searcher 


will  jump  towards  the  peak  of  the  Gaussian.  Increasing 
the  computational  domain  size  will  only  serve  to  worsen 
the  effect  shown  in  Figure  4  (right);  the  searcher  can 
wander  very  far  from  the  center  of  the  target  distribu¬ 
tion  before  returning  to  the  concentrated  area  of  targets. 

We  conclude  that  when  searching  with  apriori  infor¬ 
mation,  biasing  the  searching  direction  has  the  most  posi¬ 
tive  effect  of  the  methods  considered.  Our  results  concern 
nondestructive  search,  which  is  the  predominant  type 
seen  in  nature  [11].  As  a  point  for  further  investigation, 
we  propose  a  simple  way  to  extend  our  algorithms  to 
destructive  searching.  Start  by  performing  a  nondestruc¬ 
tive  search  with  a  priori  target  information  in  the  form  of 
a  probability  density  P,  as  above.  As  targets  are  located 
and  destroyed,  dynamically  update  P  to  include  infor¬ 
mation  about  areas  that  have  been  searched  and  thus 
are  now  known  to  contain  no  targets.  In  this  way  a  bi¬ 
ased  search  would  continue  to  look  toward  new  areas  that 
have  not  been  searched.  The  conclusions  about  bias  be¬ 
ing  more  effective  than  the  dynamic  a  search  are  largely 
based  on  our  computational  results  for  single  point  Gaus¬ 
sian  probability  distribution.  As  another  point  of  further 
study,  it  would  be  interesting  to  see  how  the  different 
methods  perform  for  bimodal  and  more  complex  distri¬ 
butions  of  target  sites.  Moreover,  a  recent  study  [12] 
considers  the  dynamics  of  Levy  flights  in  the  presence  of 
external  deterministic  potentials.  It  would  be  interesting 
to  consider  optimal  search  strategies  given  an  external 
field  or  flow.  This  is  particularly  relevant  to  underwa¬ 
ter  searching  which  is  affected  by  ocean  currents  and  to 
aerial  searching  affected  by  wind  and  weather  patterns. 
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